
%本代码的主要作用是想让训练数据的输出是5*5数据，且数据都是一样的，从而构造网络使用
%但是效果不理想

clc
clear
load 413559000.mat

% AIS.MMSI=[];
% AIS.Status=[];
% AIS.Destination=[];
% AIS.Draught=[];
% AIS.ShipName=[];
% AIS.ShipType=[];
% AIS.Length=[];
% AIS.Breadth=[];
% AIS.Time=[];
% ship_AIS = AIS;
% AIS.time=[];
% AIS.ID=[];

%速度小于1时，默认为锚泊状态，不处理
index=AIS.Speed<1.0;
AIS(index,:)=[];
count_col=size(AIS,2);
AIS(130:end,:)=[];
% 创建图形，但不显示在屏幕上
figure(1);
DrawENC(AIS, 'Original AIS Trajectory');
hold on
% plot(error_position_ais.Lon,error_position_ais.Lat,'ro');
%% 数据集划分
%将原始数据直接划分为训练集和测试集，然后在分别对其进行重构
count_ais=size(AIS,1);
idx_train=1:floor(0.9*count_ais);
idx_test=floor(0.9*count_ais)+1:count_ais;
data_train=table2array( AIS(idx_train,:));
data_test=table2array(AIS(idx_test,:));
count_row_train=size(data_train,1);%训练集的行数
count_row_test=size(data_test,1);%测试集的行数
%分别重构训练集和测试集，时间步长为5
time_step=5;
train_x=cell(count_row_train-time_step,1);
train_y=train_x;%zeros(size(train_x,1),count_col);

test_x=cell(count_row_test-time_step,1);
test_y=zeros(size(test_x,1),count_col);

%重构训练集
for i=1:count_row_train-time_step
    train_x{i}=data_train(i:i+time_step-1,:);
    train_y{i}=[data_train(i+time_step,:);data_train(i+time_step,:);data_train(i+time_step,:);data_train(i+time_step,:);data_train(i+time_step,:)];
end


%重构测试集
for i=1:count_row_test-time_step
    test_x{i}=data_test(i:i+time_step-1,:);
    test_y(i,:)=data_test(i+time_step,:);
end

%% 数据归一化,是每个维度数据的归一化
muX = mean(cell2mat(train_x));
sigmaX = std(cell2mat(train_x),0);

muY = mean(cell2mat(train_y));
sigmaY = std(cell2mat(train_y),0);
%训练数据归一化
for n = 1:numel(train_x)
    train_x{n} = (train_x{n} - muX) ./ sigmaX;
    train_y{n} = (train_y{n} - muY) ./ sigmaY;
end

%测试数据归一化
for n=1:numel(test_x)
    test_x{n} = (test_x{n} - muX) ./ sigmaX;
    % test_y(n,:) = (test_y(n,:) - muY) ./ sigmaY;
end

%% 定义LSTM网络结构
num_Features =size(train_x{1},2);   %输入特征维度
num_Responses = 5;%size(train_y,2);  %输出特征维度   可以对比单维度输出和5个维度一起输出时的结果
layers = [
    sequenceInputLayer(num_Features)
    lstmLayer(128,OutputMode="last")
    fullyConnectedLayer(num_Responses)];

%% 指定训练选项
options = trainingOptions("adam", ...
    MaxEpochs=100, ...
    SequencePaddingDirection="left", ...
    Shuffle="every-epoch", ...
    Plots="training-progress", ...
    Metrics="rmse",...
    Verbose=false);
%TargetDataFormats="CB", ... % 根据目标数据的实际格式设置
%% 训练网络
% y1=cellfun(@(x) x(:, 1:5), train_y, 'UniformOutput', false);
% y1=cell2mat(y1);
% y2=cellfun(@(x) x(:, 2), train_y, 'UniformOutput', false);
% y2=cell2mat(y2);
% y=[y1 y2];


net = trainnet(train_x,train_y,layers,"mse",options);
    %% 测试网络
    YTest = minibatchpredict(net,test_x, ...
        SequencePaddingDirection="left", ...
        UniformOutput=false);


% for i=1:2  %分别预测经纬度数值
%     net = trainnet(train_x,train_y(:,i),layers,"mse",options);
%     %% 测试网络
%     YTest = minibatchpredict(net,test_x, ...
%         SequencePaddingDirection="left", ...
%         UniformOutput=false);
%     data_predict(:,i)= cell2mat(YTest);
%     net = resetState(net);
% 
% end

%% 将获得的预测经纬度进行还原
data_predict=data_predict.*sigmaX(1:2)+muX(1:2);

plot(data_predict(:,1),data_predict(:,2),'k-o');
hold on
plot(test_y(:,1),test_y(:,2),'r*');

% % 在直方图中可视化均方误差。
% figure
% TTest = cellfun(@(x) x(:, 1), test_y, 'UniformOutput', false);
% TTest=cell2mat(TTest);
% YTest=cell2mat(YTest);
% histogram(mean((TTest -YTest ).^2,2))
% xlabel("Error")
% ylabel("Frequency")
%
% % 计算总体均方根误差。
% rmse = rmse(YTest,TTest);
%
% % 绘制预测频率对实际频率的图。
% figure
% scatter(YTest,TTest, "b+");
% xlabel("Predicted Frequency")
% ylabel("Actual Frequency")
% hold on
%
% m = min(y1);
% M=max(y1);
% xlim([m M])
% ylim([m M])
% plot([m M], [m M], "r--")